function mat = fn_expmx(nv,np)

% nv = number of variables, no limit
% np = degree of polynomials {1,2,3}

% constant and 1st-order
num = 1 + nv;

% 2nd-order terms
if np > 1
    num = num + nchoosek(nv,1); % one variable squared
    if nv > 1
        num = num + nchoosek(nv,2); % two variables
    end
end

% 3rd-order terms
if np > 2
    
    num = num + nchoosek(nv,1); % one variable to the third
    if nv > 1
        num = num + nchoosek(nv,2); % one variable squared, another variable
        num = num + nchoosek(nv,2); % one variable, another variable squared
    end
    if nv > 2
        num = num + nchoosek(nv,3); % three variables
    end

end

mat = zeros(num,nv);

row = 1;

%% put exponents
    
% 1st-order terms
for v=1:nv
    row = row+1;
    mat(row,v) = 1;
end
    
% 2nd-order terms
if np > 1
    for v1=1:nv
        for v2=v1:nv
            row = row+1;
            mat(row,v1) = mat(row,v1) + 1;
            mat(row,v2) = mat(row,v2) + 1;
        end
    end            
end
    
% 3rd-order terms
if np > 2
    for v1=1:nv
        for v2=v1:nv
            for v3=v2:nv
                row = row+1;
                mat(row,v1) = mat(row,v1) + 1;
                mat(row,v2) = mat(row,v2) + 1;
                mat(row,v3) = mat(row,v3) + 1;
            end
        end
    end            
end

return